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Abstract 

In many applications of frequency estimation, the frequencies of the signals 
are so high that the data sampled at Nyquist rate are hard to acquire due to 
hardware limitation. In this paper, we propose a novel method based on sub¬ 
space techniques to estimate the frequencies by using two sub-Nyquist sample 
sequences, provided that the two under-sampled ratios are relatively prime in¬ 
tegers. We analyze the impact of under-sampling and expand the estimated 
frequencies which suffer from aliasing. Through jointing the results estimated 
from these two sequences, the frequencies approximate to the frequency compo¬ 
nents really contained in the signals are screened. The method requires a small 
quantity of hardware and calculation. Numerical results show that this method 
is valid and accurate at quite low sampling rates. 
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1. Introduction 

Frequency estimation has wide applications in communications, audio, medi¬ 
cal instrumentation and electric systems iii- Frequency estimation methods 
cover classical modified DFT [^, subspace techniques such as MUSIC 0 and 
ESPRIT 0 and other advanced spectral estimation approaches [^. In general, 
the sampling rate of the signal is required to be higher than twice the highest 
frequency (i.e. Nyquist rate). However, it is challenging to build sampling hard¬ 
ware when signal bandwidth is large. When the signal is sampled at sub-Nyquist 
rate, it leads to aliasing and attendant problem of frequency ambiguity. 

A number of methods have been proposed to estimate frequency with sub- 
Nyquist sampling. To avoid the frequency ambiguity, Zoltowski proposed a 
time delay method which requires the time delay difference of the two sampling 
channels less than or equal to the Nyquist sampling interval Q . By introducing 
properly chosen delay lines, and by using sparse linear prediction, the method in 
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@ provided unambiguous frequency estimates using low A/D conversion rates. 
The authors of 0 made use of Chinese Remainder Theorem (CRT) to over¬ 
come the ambiguity problem, multiple signal frequencies often need the parame¬ 
ter pairing. Bourdoux used the non-uniform sampling to estimate the frequency 


ll| . Some scholars used multi-coset sub-Nyquist sam plin g with different sam¬ 


pling rates to obtain unique signal reconstruction 01 0 . Based on emerging 


compressed sensing theory, sub-Nyquist wideband sensing algorithms and corre¬ 
sponding hardware were designed to estimate the power spectrum of a wideband 


signal 0 0 16j. However, these methods usually require much hardware and 


complicated calculations, which makes the practicability discounted. 

In this paper, we propose a new method to estimate the frequency compo¬ 
nents from under-sampled measurements. The method uses two sets of sub- 
Nyquist samples and requires less hardware than former methods. When the 
sampling rates satisfy certain conditions, the frequencies can be estimated ac¬ 
curately. The paper is organized as follows: Section II describes the formulation 
of the problem and reviews ESPRIT. Section III gives our method and analysis. 
Simulation results are shown in Section IV. The last section draws conclusions. 


2. Preliminaries 

2.1. Problem Formulation 

Consider a signal containing K frequency components with unknown con¬ 
stant amplitudes and phases, and additive noise that is assumed to be a zero- 
mean stationary complex white Gaussian random process. The samples of the 
signal at sampling rate /s = 1 /At can be written as 

K 

x{n) = ^ + e(n), n = 1, 2, • • • , (1) 

fc=i 

where fk is the k-th frequency, Sk is corresponding complex amplitude, and 
e(n) is additive Gaussian noise. Assume that the upper limit of the frequencies 
is known, but we only have low-rate analog-to-digital converters whose sam¬ 
pling rates are far lower than Nyquist rate. Many articles use multi-channel 
measurement systems to estimate the frequencies. We shall demonstrate that 
double-channel sub-Nyquist sampling with specific rates is enough in general. 

2.2. ESPRIT 

The ESPRIT utilizes the subspace approach and estimates signal parameters 
via rotational invariance techniques 0. Define y{n) = x(n + 1) and choose an 
m > K. Two TO X 1 vectors x{n) and y{n) are defined as 

x{n) = [x{n), • • • , x{n + m — 1)]^, 
y{n) = [x{n -I-1), • • • ,x(n -|- to)]”^. 
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By using (CD, the time series x(n) and y(n) can be written in matrix notation, 
namely 


x{n) = As + e{n), 
y{n) = + e(n + 1), 

where s = • • • , ^ is a if x 1 vector of complex am¬ 

plitudes, e{n) = [e(n), • • • , e(n -\-m — 1)]^, and is a diagonal K x K matrix 
containing the relative phase between adjacent time samples for each of the K 
components 

$ = diag • • • , . (4) 

$ is a diagonal matrix relating the temporally displaced vectors x and y and 
is therefore referred to as a rotation operator. A is the m x K Vandermonde 
matrix, namely 


A = 


1 


1 

gj27r/KAt 


gi2ir/i(m-l)At ... gj27r/K (m-1) Ai 


The auto-covariance matrix of x is given by 

R,, = E {x{n)x{nf^ = ASA^ + a^J, 


( 5 ) 


( 6 ) 


where tr^ is the variance of e(n), I is an identity matrix and S' is a nonsingular 
matrix. (*)^ denotes the Hermite transpose operation. The cross-covariance 
matrix of the vectors x and y is given by 


Rxy = E {x{n)y{n)^^ 


=AS^^A 


H aH 


(P-Z, 


( 7 ) 


where Z is a, K x K matrix with ones on the first sub-diagonal and zeros 
elsewhere, i.e. 

“0 0 ••• 0 

1 0 ••• 0 


Z = 

0 0 10 

Rxx has (to — K) minimum eigenvalues all equal to cr^. Define 


( 8 ) 


Cxx = Rxx - cr'^I = ASA^, 

Cxy = Rxy - = AS^^A^, 

and consider the matrix pencil given by 

Cxx - iCxy = AS {I -A^. (10) 
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By inspection, the column space of ASA^ and AS^^ A^ is identical. If 
y = k-th. TOW of {I — 7 $^) will be zero and the pencil {Cxx — 

iCxy) will decrease in rank. Therefore, = 1, • • • ,K) are exactly the 

generalized eigenvalues of the matrix pair {Cxx,Cxy\- Once the generalized 
eigenvalues 7 are calculated, the frequencies of the signal can be estimated. 
However, when the sampling rate is lower than Nyquist rate, only a series of 
probable frequencies are obtained via ESPRIT. In next section, we analyze the 
problem and determine the real frequencies by sampling at two specific sub- 
Nyquist rates. 

3. Proposed Method 

3.1. The Impact of Under-sampling 

Under-sampling usually leads to aliasing of frequency spectrum. Suppose 
the highest frequency component in the signal is lower than fn, we sample at 
the rate fs = Ih/p {p> !)■ So the samples can be written as 

K 

^ -I- e(n), n = 1, 2, • • • . (II) 

fe=i 

In the case of under-sampling, if we deal with the samples as the same with that 
in normal sampling, the normalized frequency p ■ fk/fn may be greater than I 
for some k. Due to the periodicity of trigonometric functions, one 
will correspond to multiple fk. In the following we shall expand the frequencies 
within fn. 

3.2. The Expansion In Frequency Domain 

As a super-resolution subspace technique, ESPRIT has favourable robustness 
to noise. So we choose ESPRIT as the first step of our method. Then the 
generalized eigenvalues 7 in (fTUl) turn into 7 = . Beca.use p-fk/fn > I 

for some k, a series of estimated frequencies fk all satisfy = 7 . 

We restrict the under-sampled ratio p to be a positive integer. Denoting the 
principal argument of 7 as Arp ( 7 ), we obtain the estimation from sub-Nyquist 
samples 

fk = Arg{j) ■ /ij/(27rp), fc = I, 2, • • • ,K. (12) 

The fk are only a part of possible frequencies, which occupy the bottom in 
frequency domain. All the possible frequencies (we call eligible frequencies) can 
be expanded as 

fk = fk + lk- fn/p, Zfe = 0,1, • • • ,p - 1. (13) 

The eligible frequencies include the frequencies approximate to the frequency 
components really contained in the signals (we call real frequencies) and the false 
frequencies caused by under-sampling. It is worth mentioning that different 
real frequencies may correspond to the same fk- This condition has no crucial 
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influence on the following algorithms because no real frequencies will be left out 
from the eligible frequencies. Since this condition is not common, we assume 
that different real frequencies correspond to different fk for the concision of 
discussions. 

3.3. The Screening of The Eligible Frequencies 

In order to determine the true value of Ik for each k, we sample the other 
time series at another different rate f'g = fn/q, where q is another positive 
integer. How to screen the real frequencies from the eligible frequencies? An 
intuitional method is using ESPRIT for the other sample sequence again and 
then matching the two series of eligible frequencies. When the number of the 
eligible frequencies is not large, we can match the two sets of eligible frequencies 
by trial and error. That is to say, we may solve proper Ik and from a dualistic 
linear indeterminate equation 

fk + lk- Jh/p = /fc + I'k ■ fn/q, (14) 

i.e., 

qik - pi'k = pq (fk - fk) /fn, (15) 

where Ik = 0 , 1 , • ■ • — 1 and = 0 , 1 , • ■ • , g — 1, (*)' denotes the parameter 

of the second sampling. According to Bezout’s identity, when p is coprime to q 
(i.e. p J- q), (fT^ has integer solutions as long as the right hand side of (ITSl) is 

an integer. Moreover, since h G {0,1, • • • ,p — 1} and G {0,1, • ■ • , g — 1}, the 

equation (11511 just has an unique satisfactory solution. Denoting the true values 
of Ik, I'k by Ik, I'k, according to (fT^ we have 

pq (fL - fn) /fn = pq {fm - fn+ L ' fn/p - C ' fn/q) /fn 
= pq {fm - fn)/fn + qin “ pl'm- 

To make the matching correct, m must be an integer when and only when 
m = n. When in ^ n and A/ = fm — fn is the integer multiples of fn/ ipq), 
the matching of the frequencies may be suffer from mistakes. Fortunately, in 
the simulations we find that the probability of this situation is quite low. The 
condition that under-sampled ratios p and q are coprime offers good properties 
for screening, as we will see later. We can match the eligible frequencies by 
graphic method. The eligible frequency points fk and /(, are plotted together 
and then the coincident points are found as the real frequencies. When there 
are too many eligible frequencies, the direct matching of the eligible frequencies 
is impractical. In the following we propose a more practical algorithm to screen 
the real frequencies. 

We use a MUSIC-like algorithm to determine the real frequencies instead of 
the matching process. The second sample sequence can be written as 

K 

x'{n) = -b e'(n). (17) 
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Let gk = fk/fn, the time series of second sampling is expressed as 

x'{n) = [a{giJi), a{g 2 , h), ■■■ , a{gK, Ik)] s + e'(n), (18) 

where 

^{gkj ^k) — 

and m' is the number of sequential samples. The subspace spanned by a{gk, Ik) 
can be expanded. That is to say, for each k, a{gk, Ik) are augmented as 


J'^'^q(gk+lk/p) ... pj'^-^qrn'(gk+lk/p) 


(19) 


^(y9ki ^k) — ^oifk)i ^l(/fc), * ‘ * 5 l(/fc) 


( 20 ) 


where 

^rigk) — 

So (IT^ changes into 


g^32Trq(gk+rlp) . . . (gk+rjp) 


x'(n) = A' s' + e'{n), 


( 21 ) 

( 22 ) 


where 


and 


A' = [ao{gi),--- ,ap_i{gi),-■ ■ ,ao{gK),--- ,ap-i{9K)], (23) 


(24) 

The following procedure is similar to MUSIC. Taking the eigen-decomposition 
of the auto-covariance matrix of x' renders 


- [C^s, Llg 


'A O' 


r 1 

o 

to 


\ 


(25) 


where the column vectors of Us and Ue are, respectively, the eigenvectors that 
span the signal subspace and the noise subspace of Rx'x' with the associated 
eigenvalues on the diagonals of A and So the pseudo-spectrum with respect 
to fk = gk ■ Jh + Ik ■ Ih/p can be given by 


PMuifk) 




{gk)UeU^ar{gk) 


fc = l,2, 


,K,r = 0,1,■■■ ,P-1- 


(26) 


The maximum K frequencies corresponding to the highest K wave peaks or the 
frequencies that have the power larger than a threshold are the final estimated 
frequencies fk- This algorithm has less amount of computations than conven¬ 
tional MUSIC because the number of eligible frequencies is generally less than 
that of uniform grid points. In (l23l) . if m' is an integral multiple of p, the mode 
vectors a^-^gk) corresponding to the same fk are mutually orthogonal, which 
makes it easy to distinguish the real frequencies. So the condition that p and q 
are required to be relatively prime provides high discriminability for screening. 
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3.4- The Procedure of The Method 

We summarize the procedure of proposed method as shown in Fig.[T] Firstly 
one sample sequence is used to obtain the eligible frequencies via ESPRIT. 
Then the other sample sequence is used to screen the eligible frequencies via the 
MUSIC-like algorithm. The maximum K wave peaks indicate the real frequen¬ 
cies, or a proper threshold is set to filter the false frequencies. 



Figure 1: The procedure diagram of the method. 


4. Simulation Results 

In this section, we shall create satisfactory signals and estimate the frequency 
components by the proposed method. The signals are composed of K complex¬ 
valued sinusoids. The frequencies are assumed to be lower than IQQHz and 
the minimum interval is larger than O.lHz. The phase angles corresponding to 
the frequency components are set to be randomly distributed in [0, 27r). In the 
simulations, the two sub-Nyquist sequences are sampled at the rate 100/5Rz 
and lOO/TRz (i.e. p = 5,q = 7). 

In the first simulation, we demonstrate the process of our method in noiseless 
condition. The number of frequency components is set to be iF = 2 and the 
amplitudes are 0.8. As shown in Fig. [5J we plot the true frequencies /fc, the 
estimated frequencies fk via ESPRIT, the eligible frequencies ft and the pseudo¬ 
spectrum of the eligible frequencies. The true frequencies are 25Hz and 50Hz. 
The number of eligible frequencies are p ■ K = 10. We can see that the eligible 
frequencies all differ by Ik ■ fn/p from the true frequencies. 

In the second simulation, we set AT = 3 and the frequencies are set to be 
lOHz, 25Hz and 50Hz respectively. The white Gaussian noise at SNR=10dB is 
added. We use the first step of ESPRIT for the sequence at sampling rate 100/p 
and the second step of the MUSIC-like algorithm for the sequence at sampling 
rate 100/q to draw one pseudo-spectrum. Then we reverse the two sequences 
to draw the other. As shown in Fig. [31 the two pseudo-spectrums indicate the 
same result. 

Then we compare the accuracy of the method with conventional ESPRIT. 
The same signals are sampled at lOOHz for the conventional ESPRIT. The 
number of the components is set to be 3 and the SNR varies from OdB to 30dB. 
The amplitudes are set to be randomly distributed in [0.1,1]. The mean square 
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(a) 


(b) 




(c) (d) 

Figure 2: The comparison of the frequencies, (a) (b) (c) f^; (d) f^. 
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Figure 3: The pseudo-spectrum with respect to fk via the MUSIC-like algorithm. 




























error of estimated frequencies different from true frequencies is computed by 


MSE = 


K 


E (/fc - fkY 


K. Average values are obtained by 500 experiments 


k=l 


for each SNR. As shown in Fig. |4l the proposed method mainly has the same 
accuracy with conventional ESPRIT, except when SNR=0dB. 



SNR 


Figure 4: The comparison of errors at different SNR levels. 


5. Conclusion 

In this paper, we propose a sub-Nyquist method to estimate the frequencies 
of frequency-sparse signals. Two sample sequences are sampled at two spe¬ 
cific rates. Through the ESPRIT algorithm for one sequence and the screening 
algorithm for the other, the real frequencies are determined. The MUSIC-like al¬ 
gorithm is proven to be valid by numerical experiments. The proposed method 
requires less hardware than other methods and holds enough accuracy. This 
method can be used in many applications such as harmonic analysis, direction- 
of-arrival estimation, etc. 
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